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ABSTRACT 


For  accnrate  frequency  estimation  Principal  Component  Autoregressive  (PC-AR) 
spectral  estimation  methods  have  received  considerable  attention  in  the 
recent  literature.  Explicit  coarautation  of  the  Eigen-decomposition  of  the 
autocorrelation  matrix  is  required  to  obtain  the  PC-AR  solution.  An 
alternative  approach  called  the  ‘‘Eigenvalue  filtering  method*  (EFM)  where 
the  eigenspaca  need  not  be  computed,  is  proposed  in  this  paper.  The  proposed 
method  utilizes  the  geometry  of  the  distribution  of  the  eigenvalues  in  a 
matrix  function  so  that  it  closely  approximates  the  pseudoinverse  of  the 
autocorrelation  matrix.  It  is  shown  via  computer  simulation  that  compared 
with  the  Forward/Backward  method,  the  proposed  method  enhances  the 
threshold  in  SNR  by  abont  6-8  dB.  Further  improvement  is  obtained  by  a 
simple  subset  selection  method  and  a  second  eigenvalue  filtering  iteration. 
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inisopuction 

Estimation  of  parameters  of  multiple  sinusoidal  signals  imbedded  in 
white  noise  is  a  very  well  researched  problem  in  signal  processing.  It  is 
also  one  of  the  fundamental  problems  encountered  in  a  variety  of  seemingly 
unrelated  fields,  e.g,  geophysics,  economics,  radar,  astronomy,  sonar,  etc.. 
Ideas  and  techniques  emerging  from  different  fields  have  transgressed  and 
intermingled,  resulting  in  a  certain  maturity  of  this  problem.  Many 
researchers  have  reported  various  aspects  of  this  problem  in  recent  and 
earlier  literature.  Useful  references  are  listed  in  Eumaresan’s  thesis  [1] 
and  also  in  a  recent  book  by  Haykin  et  al  [2]  and  in  the  forthcoming  book  by 
Kay  [31 . 

One  of  the  ideas  which  has  received  considerable  attention  is  the  use 
of  Singular  Value  Decomposition  <SVD)  of  a  large  order  (larger  than  number 
of  sinusoids)  autocorrelation  matrix  estimated  from  data  or  equivalently, 

Grt 

the  principal  component  analysis  of  the  autocorrelation  matrix.  The  SVD 
approach  produces  orthogonal  eigenspaces  consisting  of  the  signal  subspace 
and  the  noise  snbspace.  The  many  variants  of  SVD  based  frequency  estimation 
methods  that  are  presently  available,  usually  make  use  of  either  the  signal 
snbsnace  or  the  noise  snbspace  to  estimate  the  frequency  components.  One 
disadvantage  of  these  approaches  is  thb  burden  of  computing  the  eigenspaces 
explicitly.  To  avoid  this,  a  different  approach  has  been  undertaken  in  the 
present  work,  where  the  knowledge  of  the  geometry  of  the  distribution  of 
eigenvalues  has  been  exploited.  An  appropriate  function  of  the 
autocorrelation  matrix  is  suggested  to  preserve  the  signal  snbspace  and 
nullify  the  noise  snbspace.  It  is  shown  that  if  the  geometry  of  the 


distribution  of  tbe  eigenvalues  is  exactly  known,  then  tbe  suggested  matrix 
function  eliminates  tbe  noise  subspace  while  keeping  the  signal  eigenvalnes 
intact.  In  oractical  situations,  however,  the  exact  distribution  of  the 
eigenvalnes  will  not  be  known.  A  few  theorems  which  relate  the  elements  of 
a  hermit ian  matrix  to  its  eigenvalnes  have  been  invoked  to  obtain  bounds  on 
the  signal  eigenvalnes.  Then  the  estimated  bound  (or  threshold)  is  utilised 
in  the  proposed  awtrix  function  to  remove  the  effects  of  most  of  the  noise 
eigenvalnes.  Siamlations  on  the  same  data  with  a  large  number  of 
independent  noise  realizations  have  been  performed  and  the  results  have  been 
coaqsared  to  those  obtained  by  using  the  modified  covariance  or 
Forward/Backward  (FB)  method  and  the  Craaer-Rao  (CR)  lower  bound.  Without 
computing  the  eigenspace  explicitly,  the  proposed  method  enhances  the 
threshold  (of  Signal  to  Noise  Ratio)  by  about  (-8  dB  over  the  FB  method. 
Then  a  subset  selection  method  is  used  which  improves  the  threshold  by  2-6 
dB  for  the  FB  and  the  proposed  methods.  For  the  present  method  a  second 
iteration  is  then  shown  to  enhance  the  threshold  further. 

This  paper  is  arranged  as  follows.  In  Section  I  the  problem  is 
formulated  for  multiple  sinusoids  in  noise  and  a  brief  discussion  of  the  FB 
method  and  the  Principal  Component  Autoregressive  (PC-AR)  method  is 
presented.  In  Section  II,  the  proposed  eigenvalne  filtering  method  is 
described  with  theoretical  analysis  and  simulation  results.  Section  III 
consists  of  some  concluding  observations. 


SECTION  I  :  PROBLEM  FORMULATION 


la  :  Problea  Definition 

Gives  N  tansies  of  observation  data  composed  of  s  complex  sinusoids  in 
white  noise  denoted  by 

Z  j(2«f  n  +  f  ) 

x(n)  »  >  A.  e  +  x(n)  (1) 

^  n-  0,1.... N-l 

where  k-^.k^...,^  are  the  amplitudes  are  the  frequencies  and 

are  the  phases  associated  with  each  sinusoid,  the  problem  is  to 
estimate  A1,A2....Ap,  f^fj....^  and  fx.f The  observation  noise 
z(n)  is  assumed  to  be  zero  mean  and  white.  Usually,  the  number  of  complex 
sinusoids  p  is  known  or  an  estimate  of  it  is  available  and  the  phases 
»r®  assumed  to  be  random  (or  fixed  but  unknown)  and  uniformly 
distributed  between  0  and  In. 

Ib  :  Previous  Results 


The  oldest  and  probably  the  most  widely  nsed  frequency  estimator  is  the 
oeriodogram.  It  can  be  shown  [8]  that  for  u  *  1,  periodogram  maximizes  the 
likelihood  function.  But  for  p  >  1,  the  interaction  between  the  frequencies 
causes  poorer  estimates  except  when  the  frequencies  are  far  enough  apart  or 
if  they  are  separated  by  C/N  ,  where  K  is  an  integer  not  equal  to  a  multiple 
of  N.  Especially,  for  closely  spaced  sinusoids,  i.e,  when  the  separation 
between  two  sinusoids  is  less  than  1/N,  then  oeriodogram  cannot  resolve  the 
two  peaks  and  exhibits  a  single  peak  instead  of  two.  The  maximum  likelihood 
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estimator  for  p  greater  than  1  is  ooarautationally  intensive  whether  grid 
seareh  [1]  or  iterative  techniques  [3,4]  are  nsed.  High  resolution  spectral 
estimation  techniques  have  been  extensively  nsed  by  many  researchers  to 
overcome  the  problems  of  the  periodogram.  Since  the  approach  taken  in  this 
paper  is  an  indirect  implementation  of  the  principal  component 
autoregressive  frequency  estimator  [S],  a  brief  review  of  this  method  is  in 
order.  The  proposed  method  is  then  described  in  Section  II. 

The  data  composed  of  p  sinusoids  as  given  by  equation  (1)  may  be 
thought  of  as  a  limiting  form  AR  process  with  narrow  band  peaks  at  the  p 
sinusoidal  frequency  locations.  Such  an  AR  process  can  be  modeled  as  the 
output  of  an  AR  filter  driven  by  white  noise.  With  this  model  at  hand,  one 
can  use  a  reasonable  estimator  of  the  autocorrelation  function  (ACF)  from 
the  data  and  then  estimate  the  AR  parameters  from  the  following 
relationship, 

a  -  -  R"1  r  (2) 


The  last  equation  is  the  solution  of  the  following  normal  equations 
encountered  in  FB  method. 


(3) 


or  R  a  *  -  r 


Safer  >THfe 


where , 


1  r  *■!  N-L  1 

'II  *  -  I  ^  x(n-j)  +5  x(n+i)  x(n+j)  I,  1  _  J’  *  ’  '  » 

iJ  2(N-L)  I  “  \  J  j  o,  l.  .  .  .  L. 


•re  the  ACF  estimates  used  io  the  FB  method.  The  frequency  estimates  ere 


obtained  hy  choosing  p  zeros  of  A(z)  that  are  closest  to  the  unit  circle. 


where  A(x)  *  1  +  £  aj  z~*  is  known  as  the  estimated  prediction  error  filter 


(PEP).  A  Choice  of  L  snch  that  N/3  <  p  <  N/2  has  been  shown  to  yield 


reasonably  Rood  estimates.  At  low  SNR  valnes  (below  20  dB).  the  FB  method 


fails  to  resolve  closely  spaced  peaks  dne  to  the  effect  of  noise  on  the  ACF 


estimates.  This  problem  has  been  alleviated  to  a  great  extent  by  the  use  of 


the  principal  component  solution.  In  the  PC-AR  method,  the  following  eigen- 


decomposition  of  the  R  matrix  is  used. 


1  xi  i"  *  1 


i~o+i 


where,  Xj  2  kj  -  '  ‘  -  *p+l  -  *  *  -  are  ****  eigenvalnes  (real  valued 


and  positive,  because  R  is  hermitian  and  positive  definite)  and  v^,  vj. 


▼L  are  the  corresponding  orthonormal  eigenvectors.  Using  this  notation,  the 


solution  in  (2)  can  now  be  expressed  as. 


P  L 

T  1  H  V  1  H 

£  \ ,-i  Zi  1  l  X  Zi  Zi  1 

i-1  i=p+l  1 


A  principal  component  solution  would  omit  the  second  term  which  corresponds 


to  the  L— o  noise  eigenvalues  and  retain  the  first  term  that  corresponds  to 


the  p  signal  eigenvalnes  so  that  the  PC  solution  is  given  by. 
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5PC 


-1 


1  H 

-i  -i  - 


(5) 


Equation  (5)  it  also  vrittan  as. 


-PC 


-  E  r 


(6) 


where  R^,  which  is  commonly  known  as  the  pseudo  inverse  of  R,  is  defined  as. 


_#  Y  1  H 

*  ‘/T!ili 


(7) 


i=l 


As  an  interesting  interpretation  of  the  ap£  solution  that  should  he  noted 

is  that  if  R  is  singular  the  normal  equations  (3)  will  have  infinitely  many 

solutions  out  of  which  one  unique  solution  that  minimizes  the  Euclidean  norm 
H 

a  a  is,  in  fact,  the  ap^  solution  given  by  (6).  This  is  also  known  as  a 
'miniaram-norm'  solution.  For  the  FB  case,  this  occurs  when  the  order  L  is 
greater  than  the  number  of  sinusoids  p  and  no  additive  noise  is  present. 
Also  for  the  noiseless  FB  oase,  the  corresponding  PEF  A(z)  has  p  zeros  on 


the  unit  circle  at  the  sinusoidal  frequency  locations  and  the  other  L-n 
zeros  of  the  PEF  lie  inside  the  unit  circle  [1], 


SECTION  II  :  NET  METHOD 


II*  :  Motivation 

The  PC-AR  solution  retains  the  p  smallest  eigenvalues  of  R~*  and  zeros  out 
the  L— p  largest  eigenvalues.  This  operation  can  he  functionally  represented 
as, 

*PC  =  ~  f (R_1)  r  (8) 

where  f(R~*)  defines  the  functional  operation  of  zeroing  out  the  L-p  largest 
eigenvalues  of  R_1,  i.e,  f(R_1)  -  R*,  where  R^  is  defined  in  (7).  The 
eigenvalues  corresponding  to  the  matrix  f(R-*)  are  given  by  f(A),  where  A 
denotes  the  eigenvalues  of  R-*  which  are  related  to  the  eigenvalues  of  R 
as, 

Aj  “  1/X.j  ,  A2  =  I/Aj*  •  •  •  A^  ■  1/A^  80  that  Aj  -  Aj  -  •  .  —  A^ 

In  R^  as  given  by  equation  (7),  all  eigenvalues  greater  than  A„.  i.e,  the 
L-p  largest  eigenvalues  are  zeros  which  may  be  represented  as 
f(A)  *  A  for  A  i  Ap 

-  0  for  A  >  \  (9) 

Geometrically,  f(A)  as  described  by  equation  (9),  can  be  represented  by  the 
plot  in  Fig.  1.  This  geometry  suggests  that  one  could  as  well  find  a 

function  g(A)  that  closely  approximates  f(A)  of  Fig.  1  and  then  use 

the  corresponding  matrix  function  g(R-1)  in  place  of  R^  in  the  apc 
solution  in  equation  (6).  The  advantage  of  using  g(R-*)  instead  of  R^  is 
that  the  eigen-decomposition  of  R-*  need  not  be  computed  since  g(R_1)  has 

the  same  effect  on  R~*  as  the  use  of  R^.  This  approach  whereby  a  desired 

matrix  function  f(A)  is  replaced  by  an  easier-to-calculate  matrix  function 
g(A),  where  A  is  any  L  x  L  matrix,  has  received  some  attention  in  numerical 


4 

f 
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analysis  [7].  These  techniques  are  based  on  the  idea  that  if  g(A),  the 
eigenfunction  corresponding  to  g(A).  closely  approximates  f(A)  on  the 
eigensnsce  of  the  corresuonding  matrix  A,  then  f(A)  is  annroxisiated  by  g(A). 
The  approximate  matrix  functions  do  not  involve  eigenvalues  and 
eigenvectors,  thus  avoiding  explicit  eigen-decotroosition. 

lib  :  Possible  Approximations  of  f(A) 

It  can  be  easily  shown  [7]  that  the  function  g(A)  of  an  L  x  L  matrix  A, 
has  the  same  eigenvectors  as  the  original  matrix.  Also,  if  A  is  any 
eigenvalue  of  the  original  matrix  A,  then  the  corresponding  eigenvalue  of 
the  matrix  function  g(A)  is  simply  g(A).  Some  approximating  functions  are 
now  discussed. 

1)  Polynomial  Approximation  :  The  coefficients  of  the  polynomial, 

g(A)  =  «o  +  *i  A  +  gj  A2  +  .  .  +  *m  $ 

can  be  found  by  a  least  squares  fit  with  f(A)  at  equally  (or  unequally) 
spaced  A  values.  M  may  be  any  arbitrary  integer.  This  approximation  was 
attempted  using  IMSL  Package  Subroutines  and  also  by  directly  programming  it 
in  Fortran.  But  as  shown  in  Fig. 2  for  the  case  of  M  =  10,  the  approximations 
always  yielded  a  poor  fit,  possibly  because  of  the  discontinuity  of  f(A)  at 
A  *  Aq.  Also,  g(A)  is  seen  to  have  negative  going  values  for  some  regions  of 
A  >  Ap.  The  effect  of  this  may  be  adverse  in  the  eqnivalent  matrix  function 
because  negative  values  of  g(A)  may  result  in  negative  eigenvalues  of  g(R- 
*).  In  such  cases,  the  resulting  matrix  function  will  not  be  nositive 
semidef inite,  which  is  inconsistent  with  the  property  of  a  covariance 


AV.V-V. 
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matrix.  The  frequency  estimates  obtained  using  this  aonroximat ion  technique 
were  ooor  and  hence  the  polynomial  approximation  seemed  unacceptable. 

2)  Rational  Function  Approximation  :  The  use  of  a  rational  function  of  A 
given  by. 

N(A) 

g< A)  -  -  (ll) 

D(A) 

is  now  examined.  One  choice  of  g(A)  which  closely  approximates  f(A)  in  (9) 
is  given  by, 

A 

g(A)  =  -  for  0  <  A  <  \  (12) 

1  +  (A  /  A)* 


where  A  is  chosen  such  that.  Ap  <  A  <  Ap+j.  It  can  be  easily  verified  that 
for  sufficiently  large  M 

g(A)  *  A  for  A  <  A 

»  0  for  A  >  A 

In  Fig.  3  (l/A)g(A)  vs.  A  /A  is  plotted.  It  is  evident  that  even  for  small 
values  of  N,  g(A)  dies  off  quickly  enough  for  A  >>  A.  In  practical 
situations  most  of  the  noise  eigenvalues  of  are  much  larger  than  Ap  *° 

that  even  if  a  low  value  of  M  is  used,  those  noise  eigenvalues  will  be 
substantially  suppressed. 

He  :  Analysis  of  the  rational  matrix  function  g(R~*) 


The  matrix  function  corresponding  to  the  eigenvalue  fnnction  g(A)  in  (12) 


can  be  written  es 


where  I  is  the  L  x  L  identity  matrix.  It  is  nroved  that  g(R_1)  — >  R^ 
becomes  large.  Since  the  eigenvalues  of  R~*/X  are  Aj/X.  Aj/X.  .  .  .  . 
and  R  1  and  R  1/X  have  the  same  set  of  eigenvectors  v^.vj.  ....  vL, 
*)  has  the  following  deconmosition. 


(13) 

as  N 

V* 

g(R" 


g(R  1)  •  y 

[V  «1 
• 

y® 

i  +  y 

'  -•  1 
e 

• 

y® 

• 

©i 

■ 

0  *  v 

-  < v»  J 

where  ,  V  -  [vj  v2  .  .  v^] .  Now, since  R  is  a  uositive  definite  matrix  its 
eigenvectors  form  an  orthogonal  set.  so  that,  VV®  ■  I  if  the  eigenvectors 
are  normalized  to  have  unit  lengths.  Since  R-1  has  the  same  eigenvectors  as, 
one  can  write. 


g(R_1)  -  V 


0 

V  " 


\  J 


V®V 


[l+tAj/X)11] 


0  ’  [i+o^/X)*1]  1 


v® 


A1/(l+(A1/X),,)-g(Al) 


Aj/(  l+tAj/Xj^-gtAj) 


Y/d+(\/X)M)-g(\) 


v8 
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This  derivation  verifies  the  statement  made  in  Section  lib  about  the 
eigenvalues  and  eigenvectors  of  a  function  of  a  matrix. 

Nov.  for  large  enough  M, 

g(A)  — >  0  for  A  >  X  and 
— >  A  for  A  <  X 

So  ve  have  for  large  M, 


g(R  *)  — >  V 


The  above  analysis  justifies  that  if  the  proper  scaling  parameter  X  is  known 
such  that  Ap  <  X  <  Ap+j.  then  the  suggested  matrix  function  indeed  tends  to 
for  large  M.  In  practical  situations,  however,  X  will  not  be  known  unless 
one  coamutes  the  eigenvalues  explicitly,  but  that,  in  fact,  was  intended  to 
be  avoided  in  the  first  place.  In  the  next  section,  two  matrix  theorems  are 
stated  to  obtain  a  reasonable  estimate  of  X. 


lid  :  Choice  of  X 

It  is  a  well  known  fact  that  the  trace  of  matrix  is  equal  to  the  sum  of 
the  eigenvalues  of  the  matrix.  The  noise  eigenvalues  Ap+j.  •  •  are  much 
larger  than  the  signal  eigenvalues.  Therefore,  it  can  be  expected  that  the 
average  of  the  eigenvalues  will  be  larger  than  Ap.  So  one  choice  of  X  is 
X'  *  (1/L)  •  sum  of  the  eigenvalues  of  R“*  *  (1/L)  2  R-*(i.i)  *  (1/L) 
Trace(R~*),  where  R-*(i,i)  denotes  the  (i,i)th  element  of  the  matrix  R-*. 
Also,  sinoe  R-3,  is  a  hermit isn  matrix  with  eigenvalues  Aj  1  ^  1  .  .1 


y» 
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the  following  fact  is  also  true  [6]. 

If  02.03.  .  .  ,cL  and  dj .dj.  .  .  .dL  are  real  numbers  such  that, 
d.  >  R^a+l-i.L+l-i)  +cj  I  |R-1(L+l-i,j)|2  i  =  1,  2,  .  .  ,L 
ct  >0.  d£  -  dj.j  2  l/ct  .  i  =  2,3,  .  .  .  L  (14) 

than  Aj  —  dj.  i  *  1.2,  .  .  ,L 

In  numerous  simulation  examnles  with  large  number  of  noise  realizations,  the 
explicit  computation  of  the  eigenvalues  of  R-2  indicated  that  if  Cj  *  1  for 
all  i,  i.e.  if  the  d^'s  are  chosen  such  that 

d£  »  R-1(L+l-i.L+l-i)  +  Z  |  If^L+l-i.j)  j2  i  «  1,2.  .  .  ,  L. 

j  >1-41-1 

then  Aj  1  dj  Vi.  i.o.  the  eigenvalues  were  found  to  be  always  less  than 
the  corresponding  dj's.  So  the  d^'s  obtained  from  the  above  expression  can 
be  used  as  the  upperbounds  of  the  A^'s.  For  the  simulation  examole 
described  in  the  following  section  where  p  *  2  was  considered.  d2 
provides  an  upperbound  for  A2.  So  another  choice  of  A  is.  X*  “  R  *(L-1,L-1) 
+  |R-2(L-1,L)  I2.  Finally,  X  was  chosen  as  the  smaller  of  the  above  two 
bounds,  i.e,  X*  min  (  X'  .  X*  ). 

lie;  Simulation  Results 

For  the  process  given  by  (1)  with  two  sinusoids  or 

x(n)  -  Aj  eJ(2nfln  +  *i>  +  A2  e  J(2"f2n  +  +  z(n)  (15) 

where  fj  »  0~52,  f2  ■  0^50,  -  1.0,  A2  *  1.0,  dj  «  n/4  and  “  °*  is 

a  computer  generated,  complex  white  gauss ian  noise  sequence  with  variance 
202,  where  o2  is  the  variance  of  the  real  and  imaginary  parts  of  z(n).  SNR 
for  either  sinusoid  is  defined  as  10  IoSjqIAj2^®2) .  For  every  trial  25  (=*N) 
data  points  were  generated.  This  same  data  set  was  used  in  references  [1] 


and  [51.  For  every  trial  the  sinusoidal  signal  was  kept  unchanged  and  a 
different  realisation  of  z(n)  was  generated.  500  such  trials  were  performed 
for  the  FB  method  and  the  proposed  method.  The  mean,  variance  and  MSB  were 
calculated  for  fj.  In  Table  I,  bias  and  MSE  are  tabulated  for  different  SNR 
values.  A  plot  of  lOlog^^d/MSE)  vs.  SNR  is  shown  in  Fig.  4  for  the  present 
method  (plot  marked  EFM)  and  the  FB  method.  For  the  FB  method  with  L*8  which 
exihibited  the  best  results,  the  threshold  MSE  occurs  at  about  an  SNR  of  20 
dB  and  for  the  present  method  with  L*10  and  M«6  which  produced  the  best 
results,  the  threshold  ocoured  at  14-15  dB.  As  shown  in  Fig.  4,  improved 
performance  is  obtained  for  L*10  and  M**6  with  a  second  eigenvalue  filtering 
iteration  which  will  be  explained  in  Section  Ilh.  The  Cramer-Rao  lower 
bound  for  an  unbiased  estimator  is  also  tabulated  in  Table  I  and  plotted 
in  Fig.  4  for  comparison. 

To  indicate  the  iaiorovemeut  in  performance  which  may  be  obtained  if  A  could 
be  more  accurately  determined,  the  threshold  parameter  A  was  computed 
explicitly  as  (Aq  +  Ap+j)/2.  Then,  the  present  method  exhibited  results  very 
close  to  the  PC-AR  method  when  L  *  12  to  16  was  used  with  M  =  4,5.  The  best 
results  are  obtained  for  L  *  Id  and  M  *  5  as  seen  in  Fig. 4. 

Ilf  :  Commutation  of  amplitudes  (A^'s) 

Once  the  estimates  f^,  f2  are  obtained,  the  complex  amplitudes  Acj  * 
A^e^i  and  Ac2  *  A2e^*  are  computed  as  the  least  squares  solution, 

IAci  ac2]t  -  tV  h)'1  h+  ?  <™> 

where  *T’  denotes  matrix  transpose  and  *  +  *  denotes  matrix  transpose 


conjugated  and 
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|  ^ntN-Ufj  aJ2«(N-l)f2  j 

and  j  ■  [  x(0)  x(l)  .  .  .  x(N-l)  .  The  absolute  values  of  A 

are  used  as  tbe  aaolitude  estiaates  Aj  and  A2. 


lit:  Threshold  enhaneeaent  by  subset  selection 

A  siaole  subset  selection  aethod  for  frequency  estiaates  which  performed 
well  in  siaulations  is  now  described.  Coaputationally,  this  aethod  is 
relatively  less  expensive  than  siailar  approaches  reported  earlier  in  [9] 
and  [10].  Instead  of  choosing  only  2  xeros  of  A(x)  nearest  to  the  unit 
circle.  4  zeros  closest  to  the  unit  eircle  were  chosen.  The  aaplitudes 
corresponding  to  those  4  frequency  estiaates  were  coaputed  froa. 


<  A.l  »o2  A„J  *.4  1T  -  I  ?4+  ?4  I'*  -4+  I 


where, 

B.  ■  I  e.  e  e  e.  I.  and  the  e  's  are  similarly  defined  as  in  (17). 

-4  L  -f2  -f2  -f3  -f4  J  -ri 

Then  the  frequencies  corresponding  to  the  two  largest  amplitude  estimates 
were  used  as  frequency  estiaates.  The  results  of  subset  selection  aooroach 
are  plotted  in  Fig.S.  Cornered  to  the  case  of  choosing  the  frequency 
estiaates  froa  the  two  zeros  closest  to  the  nnit  circle,  about  2-3  dB 
threshold  enhancement  is  observed.  BFM  was  used  for  L*8  and  M»15  in  both  the 
oases.  For  the  case  of  n  sinusoids,  2p  zeros  of  A(z)  closest  to  the  unit 
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circle  would  produce  2o  tmlitnda  estimates,  out  of  which  tbe  frequency 
estiastes  should  be  cbosen  corresponding  to  tbe  9  largest  aaolitudes. 
Siaulation  experience  indicates  that  instead  of  2o  zeros,  if  even  store  zeros 
olosest  to  tbe  nnit  circle  are  nicked  to  estiaate  tbe  largest  p  amplitudes, 
further  threshold  enhanceaent  aay  be  obtained.  Using  L  =  8  and  M  =  15,  tbe 
best  results  were  obtained  when  tbe  two  frequencies  bad  been  chosen  in  the 
saae  manner  froa  tbe  six  zeros  of  A(z)  closest  to  the  unit  circle.  As  shown 
in  Fig.  5,  about  3  dB  further  threshold  enhancement  is  observed. 


Ilh  :  Threshold  Enhancement  by  a  Second  iteration 

The  choice  of  A  as  described  in  section  lid  leaves  a  few  noise  eigenvalues 
unchanged  in  the  approximating  aatrix  function  g(S~*)  because  the  A  chosen 
according  to  those  bounds  is  usually  larger  than  a  few  eigenvalues  which 
are  greater  than  Ap*  I®  the  second  iteration,  one  can  eliminate  a  few  more 
noise  eigenvalues  froa  g(R-1)  following  the  saae  procedure  as  in  the  first 
iteration.  A  should  be  ohosen  as  the  average  of  the  diagonal  elements  of 
g(g-*),  because  the  first  iteration  reduces  most  of  the  noise  eigenvalues, 
causing  the  average  eigenvalue  to  be  decreased.  In  the  second  iteration,  the 
new  matrix  function  can  again  be  obtained  from  equation  (13)  where  R 
should  be  replaced  by  g(R-*)  and  the  new  A  should  be  used  for  scaling.  The 

best  results  were  obtained  for  L  -  10  and  M  *  6.  The  results  are  plotted  in 
Fig.  4  along  with  the  previous  results.  This  second  iteration  reduces  the 
MSE  so  that  the  plot  is  closer  to  the  CR  bound  line  and  the  threshold  also 
decreases  by  2-3dB.  The  subset  selection  method  suggested  in  the  section  Ilg 
provided  further  enhancement  of  this  threshold.  The  best  results  were 
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obtained  for  L  *  12  and  M  ■  4,  as  shown  in  Fig.  5.  The  olots  show  the 


results  of  two  eases.  For  the  first  ease,  the  two  frequencies  having  the  two 


largest  aaelitudes  were  chosen  out  of  the  six  zeros  of  A(z)  closest  to  the 


unit  circle,  whereas  for  the  second  case,  eight  zeros  closest  to  the  unit 


circle  were  used  to  obtain  the  two  frequency  estimates.  As  seen  in  Fig.  5., 


the  second  case  yielded  the  best  performance  compared  to  the  Cramer-Rao 


bounds  and  the  known  X  case. 


SECTION  III:  OBSERVATIONS  AND  CONCLUSIONS 


It  has  been  shown  that  if  X  is  known  exactly,  the  proposed  eigenvalue 


filtering  aMthod  works  almost  as  well  as  the  Principal  Component  based 


frequency  estimation  method.  If  X  is  estimated  from  the  covariance  matrix. 


an  improvement  of  5-6  dB  in  the  threshold  is  obtained  over  the 


Forward/Backward  method.  Additional  threshold  enhancement  of  about  6-7  dB 


can  be  obtained  with  a  second  iteration  at  the  cost  of  more  computation.  A 


PEF  order  L  «  N/3  was  observed  to  yield  the  best  results  when  only  the 


first  iteration  was  used.  A  possible  reason  for  this  may  be  that  for  L  < 


N/3,  the  signal  eigenvalues  are  too  noise  corrupted  while  for  L  >  N/3,  the 


eigenvalue  filtering  method  preserves  too  many  noise  eigenvalnes.  A  better 


choice  of  X  should  enhance  the  threshold  at  least  to  that  of  the  Principal 


Component  based  techniques.  A  simple  subset  selection  method  has  also  been 


suggested  to  taprove  the  threshold  for  the  Forward/Backward  and  the  proposed 


methods. 


Ill  ■  HI  — Bi  li  |  B  !■!  II  I  I  I  I  I 
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^9*  Geometric  representation  of  the  distribution  of  the  eigenvalues 

of  the  pseudoinverse  of  The  plot  depicts  that  the  pseudoinverse 
retains  only  the  signal  eigenvalues  At#  .  .,Ap  and  sets  the  rest 
of  the  eigenvalues  to  zero  . 


Rational  function  approximation  of  f(A).  Normalized  version  of 
g(A)  as  in  eqn.  (12)  is  shown  for  M  =  2,3,  .  .  ,15. The  desired 
f(A)  is  shewn  by  the  dotted  line.  The  rational  function  approx- 
-imation  is  seen  to  approach  the  desired  one  as  M  increases. 


Performance  carper ison  of  the  Eigenvalue  filtering 
method  with  FB  method  and  CR  bound.  Frequency  estimates 
were  chosen  from  the  two  zeros  of  the  PEF  closest  to 
the  unit  circle  on  the  z-plane.  Plots  are  shown  for  the 
estimate  of  f,(=0.52). 


o 

o 

• 

o 


SNR  (dB) 


Fig.  5.  Performance  comparison  of  the  Eigenvalue  filtering 
method  with  OR  bound.  The  two  frequencies  having  the 
largest  amplitudes  were  chosen  out  of  4-8  zeros  of 
the  PEF  closest  to  the  unit  circle  on  the  z -plane . 
Plots  are  shewn  for  the  estimate  of  f^(«0.52) . 
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